Characterization of Coastal Systems#

Welcome to the first notebook exercise of Coastal Systems (TU Delft, MSc Coastal Engineering)! This is the first year that we will experiment with notebooks in this course. With these notebooks we hope to provide you with interactive course material that helps you better understand the processes and concepts that we teach in this course. Please let us know how you find the notebooks - we appreciate your feedback!

Chapter 2 of Coastal Dynamics Open Textbook describes the large geographical variation of coasts across the world. It explains how the coasts that we have today are shaped by both present-day processes and processes millions years ago. It distinguishes between three different order of features, which are associated to different orders of of time. In this notebook we will look at coastal systems at these different orders of scale.

Import libraries that we use for our analysis#

In the two cells below we import the libraries that we need for the analysis. We also set some path settings to load the data and source code. For example, in the cell below we add the src directory to the system path, which allows us to import generic functions from ../../src/coastpy.

import os
import pathlib
import sys

# make coastpy library importable by appending root to path
cwd = pathlib.Path().resolve()
proj_dir = cwd.parent.parent  # this is the root of the CoastalCodeBook
sys.path.append(str(proj_dir / "src"))
import warnings

os.environ["USE_PYGEOS"] = "0"  # use shapely 2.0 instead of PYGEOS and silence warnings

import colorcet as cc
import dask.dataframe as dd
import geopandas as gpd
import holoviews as hv
import hvplot.pandas  # noqa: API import
import hvplot.xarray  # noqa: API import
import ipyleaflet
import numpy as np
import pandas as pd
import panel as pn
from geoviews import tile_sources as gvts
from ipyleaflet import Map, Marker, ScaleControl, basemaps
from ipywidgets import HTML

pn.extension()

from coastpy.geometries import geo_bbox

# definitly not a best practice, but a workaround to avoid many warnings triggered by shapely 2.0 release
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)

DATA_DIR = proj_dir / "data"
coastal_systems_fp = DATA_DIR / "01_coastal_systems.gpkg"

Exercise 1: Plate tectonics#

In this notebook we will start with the broadest (or first order) features of the coast that cover large geographical distances (thousands of kilometres) and are linked to the long-term geological process of plate tectonics. We will do so by using earthquake data from the USGS. The dataset we load contains a sample (10%) of observed eartquakes between Jan 2000 and Dec 2018. Why earthquake data? Earthquake data reveals geologists the mysteries of the deep, but also for coastal researchers the data is insightful. Let’s first load the data by running the next cells.

Load the earthquake data#

We load the data (tabular data including geometries) and index the columns to only keep the data in memory that we actually need. In total the dataset contains 2.2 million earthquakes, but here we use a sample (10%), so the data contains approx. 220k eartquake entries. If you find that the interactive panel responds slow to how you adjust the widgets, please consider to take another sample. You can do so by uncommenting the sample line in the next cell. So if you set frac=0.1 you have a dataframe with approx. 22k unique earthquakes over the world.

WEB_MERCATOR_LIMITS = (
    -20037508.342789244,
    20037508.342789244,
)  # max polar latitudes that can be handled in World Mercator

df = (
    dd.read_parquet(DATA_DIR / "01_earthquakes_sample.parquet")
    .sample(
        frac=0.1
    )  # uncomment this line if loading the data takes too long on your computer
    .set_index("time")
    .compute()
    .tz_localize(None)
    .sort_index()
)


# To save memory we drop most of the columns. Also we drop the polar latitudes that cannot be displayed in the web mercator projection.
df = df[["mag", "depth", "latitude", "longitude", "place", "type"]][
    df["northing"] < WEB_MERCATOR_LIMITS[1]
]
df.head()
mag depth latitude longitude place type
time
2000-01-02 18:48:55.600 3.90 32.000 35.480 22.510 central Mediterranean Sea earthquake
2000-01-03 03:43:08.900 2.05 7.399 34.725 -116.301 13km W of Ludlow, California earthquake
2000-01-03 05:47:06.280 1.40 1.602 34.272 -116.446 18km N of Yucca Valley, California earthquake
2000-01-04 08:13:20.000 2.00 11.000 38.100 -90.480 7km ESE of De Soto, Missouri earthquake
2000-01-04 20:48:47.450 1.70 9.699 34.670 -116.292 13km WSW of Ludlow, California earthquake

Visualization of the earthquake data#

To explore the data we use visualization tools from the Holoviz project that makes high-level tools to simplify visualization in Python. In the next cell we enable the interactive mode on the data dataframe, create widgets to explore the data and filter the dataframe accordingly. To explore the eartquake data we create an overlay of the eartquakes on a tileset of ESRI Imagery. Please note that the code in the next cell will only do the computations and store the result in an object called panel. To actually see the results you have to run one more cell; the one that calls this object panel.

pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Exercise 1: Plate tectonics & first-order coastal features",
        style={"color": "black"},
        width=800,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

# define widgets that can be used to index the data
magnitude_slider = pn.widgets.RangeSlider(
    name="Earthquake magnitude [Richter]", start=0.1, end=10
)
depth_slider = pn.widgets.RangeSlider(name="Earthquake depth [km]", start=0.1, end=650)
date_slider = pn.widgets.DateRangeSlider(
    name="Date", start=df.index[0], end=df.index[-1]
)
column_types = pn.widgets.Select(options=["mag", "depth"])


@pn.depends(
    magnitude_slider.param.value_start,
    magnitude_slider.param.value_end,
    depth_slider.param.value_start,
    depth_slider.param.value_end,
    date_slider.param.value_start,
    date_slider.param.value_end,
    column_types.param.value,
)
def plot_earthquake_panel(
    magnitude_start,
    magnitude_end,
    depth_start,
    depth_end,
    date_start,
    date_end,
    column_type,
):
    panel = df[(df.mag > magnitude_start) & (df.mag < magnitude_end)]

    panel = df[
        (df.mag > magnitude_start)
        & (df.mag < magnitude_end)
        & (df.depth > depth_start)
        & (df.depth < depth_end)
        & (df.index >= date_start)
        & (df.index <= date_end)
    ]
    # inverted fire colormap from colorcet
    cmap = cc.CET_L4[::-1]
    colorbar_labels = {"mag": "Magnitude [Richter]", "depth": "Earthquake depth [km]"}

    return panel.hvplot.points(
        x="longitude",
        y="latitude",
        geo=True,
        color=column_type,
        global_extent=True,
        tiles="ESRI",
        frame_width=900,
        ylabel="Latitude [deg]",
        xlabel="Longitude [deg]",
        cmap=cmap,
        tools=["tap"],
        hover_cols=["place", "time"],
        logz=True,
        clim=(1, None),
        clabel=colorbar_labels[column_type],
    )


cmap = cc.CET_L4[::-1]  # inverted fire colormap from colorcet

earthquake_panel = pn.Column(
    title_bar,
    pn.Row(column_types),
    pn.Row(pn.Column(magnitude_slider, depth_slider)),
    pn.Row(depth_slider),
    pn.Row(date_slider),
    pn.Row(plot_earthquake_panel),
)

If the visualization is too slow, please follow the instructions in loading the data for taking a sample.

After running the cell below you will have a panel with several widgets to index the eartquake data; by magnitude, depth and time, while the colors on the map show either the magintude or the depth of the earthquakes.

earthquake_panel

Explore the earthquake data & questions#

  1. How do the earthquake magnitude and earthquake depth relate to the coasts that we see? (Hint: See Figure 2.3 in the textbook and consider how deep under the ground the plates are moving. Extra hint: How do earthquake magnitude and depth differ for convergent and divergent plate boundaries?)

  2. Earthquake data support one of the most fundamental processes in the geology: plate tectonics. Although plate tectonics is a relatively slow process that acts on the geological time scale, it has had an enormous impact on the formation of coastlines and determines the broadest features of the coast. What are two important inherited aspects of this process? (Hint: see Figure 2.10 and Sec. 2.3.3 in the textbook.)

  3. In 1971 Inman, D. L. & Nordstrom, C. E. used plate tectonics to classify the coast. Explain the classification that they introduced. What are the three different classes that they distinguish? How do they match with the earthquake data as you can explore in the panel?

  4. Can you identify or predict areas around the world where you will find the coasts that are distinguished by Inman, D. L. & Nordstrom, C. E.? For instance, what kind of coasts do you have in Chili? And how are they different to the east coast of the USA? And what is characteristic about the East China sea?

  5. Inman, D. L. & Nordstrom (1971) further distinguish Afro-trailing-edge coasts and Amero-trailing-edge coasts based on differences in sediment supplies. What is the main cause of these differences in sediment supply? And how do you expect the differences in sediment input to show in the coastal geomorphology?

Exercise 2: Process-based coastal classification#

In the section part of this notebook we will explore several coastal systems around the world considering the second and third order scale of features. In chapter two of the Coastal Dynamics open textbook it is explained how coastal systems can be classified according to the processes that characterize these systems. For example, one of the figures (below) shows how the relative influence of fluvial, wave and tidal processes influences the shape of coastal features. The idea of this exercise is that you identify the signatures of the processes that are introduced in chapter 2 in several coastal systems around the world.

image

The coastal systems data#

In the cell below we define a small plot function that generates a ESRI World Imagery basemap given a longitude, latitude, zoom level and name. Also, a small sheet of coastal systems around the world is loaded into geopandas, a Python library for geospatial tabular data. In the cells afterwards we sample this dataframe and show the coastal system on a map. Since the sample is random you might encounter the same coastal system multiple times; then you can just run the cell again to get another ‘coastal draw’.

coastal_systems = gpd.read_file(coastal_systems_fp)
pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Exercise 2: Coastal system characterization",
        style={"color": "black"},
        width=800,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

options = coastal_systems.name.to_list()
coastal_systems_slider = pn.widgets.Select(
    name="Coastal system", options=options, value=np.random.choice(options)
)


@pn.depends(coastal_systems_slider.param.value)
def plot_coastal_system(name):
    system = coastal_systems.loc[coastal_systems["name"] == name].copy()
    west, south, east, north = system[
        ["west", "south", "east", "north"]
    ].values.flatten()

    return system.hvplot.points(
        x="lon",
        y="lat",
        geo=True,
        color="red",
        alpha=0,
        xlim=(west, east),
        ylim=(south, north),
        tiles="ESRI",
        frame_width=1100,
        ylabel="Latitude [deg]",
        xlabel="Longitude [deg]",
    )


app = pn.Column(
    title_bar,
    pn.Row(coastal_systems_slider),
    pn.Row(plot_coastal_system),
)
app

Explore the coastal systems#

While sampling over a range of coastal systems, try to answer the following questions.

  1. Find and compare a heavily engineered river-dominated delta and a more natural river-dominated delta

  2. Compare the scale of the biggest and smallest tidal basin in the dataset

  3. Find the estuarine and deltaic systems with a spit

  4. Compare and contrast wave-dominated deltas with high and low sediment supply. How can you tell?

  5. Find a tidal estuary with large fine (muddy) sediment supply, then find one with a large coarse (sandy) sediment supply. How can you tell the difference?

  6. The eastern and western tips of the Dutch and German Wadden Islands are very different beach ridge environments. How might differences in sediment supply explain this? Where is the sediment coming from?

  7. The Dune du Pilat in France is one of the world’s largest coastal sand dunes (it is also one of the coolest and you should definitely visit if you get the chance!). Why is it located on the east side of Arcachon Inlet and not the west?

  8. Look at the northern Jiangsu coast in China. What might explain the limited sediment supply in this location?

  9. Find an estuary or tidal bay with extensive intertidal flats. Do you see salt marshes or mangrove forests nearby? Why or why not?

  10. Find an inlet with jetties. How might this affect the way it evolves?

  11. Find a delta/estuary/inlet whose shape is constrained by the presence of rocky coastal features.

  12. The Albufeira Lagoon in Portugal opens and closes seasonally. In the image shown, is it open or closed? When and how might it open or close?

  13. Find examples of heavily urbanized estuaries. How might these human interventions influence the natural processes there?

  14. Based on these satellite images, which is the most beautiful site? Taking a moment to appreciate the beauty of these natural systems is an important part of your job as coastal engineers.

from coastpy.maps import plot_esri_basemap
m = plot_esri_basemap(-9.181333, 38.510379, 13, "Lagoa de Albufeira")
m
bbox = geo_bbox(m.west, m.south, m.east, m.north)
bbox.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
m.west, m.south, m.east, m.north
(0.0, 0.0, 0.0, 0.0)
import shapely
albufeira = {
    "name": "Lagoa de Albofeira",
    "lon": -9.181333,
    "lat": 38.510379,
    "zoom": 13,
    "west": -9.276924133300783,
    "south": 38.483560395392516,
    "east": -9.085865020751955,
    "north": 38.53729004998249,
    "geometry": shapely.geometry.point.Point([-9.181333, 38.510379])
}
albufeira = gpd.GeoDataFrame([albufeira], crs=4326)
coastal_systems = pd.concat([albufeira, coastal_systems])
coastal_systems = coastal_systems[coastal_systems["name"]!="Albufeira"]
coastal_systems.drop_duplicates().to_file(38.53729004998249)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 coastal_systems.drop_duplicates().to_file(38.53729004998249)

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/geodataframe.py:1203, in GeoDataFrame.to_file(self, filename, driver, schema, index, **kwargs)
   1120 """Write the ``GeoDataFrame`` to a file.
   1121 
   1122 By default, an ESRI shapefile is written, but any OGR data source
   (...)
   1199 >>> gdf.to_file('dataframe.shp', mode="a")  # doctest: +SKIP
   1200 """
   1201 from geopandas.io.file import _to_file
-> 1203 _to_file(self, filename, driver, schema, index, **kwargs)

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/io/file.py:535, in _to_file(df, filename, driver, schema, index, mode, crs, engine, **kwargs)
    532     df = df.reset_index(drop=False)
    534 if driver is None:
--> 535     driver = _detect_driver(filename)
    537 if driver == "ESRI Shapefile" and any([len(c) > 10 for c in df.columns.tolist()]):
    538     warnings.warn(
    539         "Column names longer than 10 characters will be truncated when saved to "
    540         "ESRI Shapefile.",
    541         stacklevel=3,
    542     )

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/io/file.py:444, in _detect_driver(path)
    442     pass
    443 try:
--> 444     return _EXTENSION_TO_DRIVER[Path(path).suffix.lower()]
    445 except KeyError:
    446     # Assume it is a shapefile folder for now. In the future,
    447     # will likely raise an exception when the expected
    448     # folder writing behavior is more clearly defined.
    449     return "ESRI Shapefile"

File ~/mambaforge/envs/codebook/lib/python3.10/pathlib.py:960, in Path.__new__(cls, *args, **kwargs)
    958 if cls is Path:
    959     cls = WindowsPath if os.name == 'nt' else PosixPath
--> 960 self = cls._from_parts(args)
    961 if not self._flavour.is_supported:
    962     raise NotImplementedError("cannot instantiate %r on your system"
    963                               % (cls.__name__,))

File ~/mambaforge/envs/codebook/lib/python3.10/pathlib.py:594, in PurePath._from_parts(cls, args)
    589 @classmethod
    590 def _from_parts(cls, args):
    591     # We need to call _parse_args on the instance, so as to get the
    592     # right flavour.
    593     self = object.__new__(cls)
--> 594     drv, root, parts = self._parse_args(args)
    595     self._drv = drv
    596     self._root = root

File ~/mambaforge/envs/codebook/lib/python3.10/pathlib.py:578, in PurePath._parse_args(cls, args)
    576     parts += a._parts
    577 else:
--> 578     a = os.fspath(a)
    579     if isinstance(a, str):
    580         # Force-cast str subclasses to str (issue #21127)
    581         parts.append(str(a))

TypeError: expected str, bytes or os.PathLike object, not float
coastal_systems2
name lon lat zoom west south east north geometry
0 Lagoa de Albofeira -9.181333 38.510379 13 -9.276924 38.483560 -9.085865 -9.085865 POINT (-9.18133 38.51038)
1 Amazon -50.158433 0.318995 9 -51.913147 -0.780005 -48.405762 1.417092 POINT (-50.15843 0.31899)
2 Columbia river -124.039699 46.247702 11 -124.477844 46.057509 -123.600998 46.437384 POINT (-124.03970 46.24770)
3 Dune de Pilat -1.217360 44.586341 13 -1.326942 44.537388 -1.107731 44.635193 POINT (-1.21736 44.58634)
4 Ebro Delta 0.869526 40.699556 12 0.650253 40.595446 1.088676 40.803675 POINT (0.86953 40.69956)
5 Ganges 90.989565 22.312137 8 87.484131 20.267350 94.498901 24.332082 POINT (90.98956 22.31214)
6 Hangzou bay 121.754241 30.488979 9 120.000916 29.537619 123.508301 31.431007 POINT (121.75424 30.48898)
7 Jiangsu coast (north) 119.829333 34.471656 12 119.610214 34.358459 120.048637 34.584888 POINT (119.82933 34.47166)
8 Jiangsu coast (south) 121.020279 32.638839 11 120.581818 32.407212 121.458664 32.869784 POINT (121.02028 32.63884)
9 Mekong river 106.762698 10.179805 11 106.324310 9.909333 107.201157 10.450000 POINT (106.76270 10.17980)
10 Missisipi -89.245992 29.151004 11 -89.684143 28.910813 -88.807297 29.390551 POINT (-89.24599 29.15100)
11 Pearl (Zhujiang River estuary), 113.636915 22.755011 12 113.417702 22.628272 113.856125 22.881553 POINT (113.63692 22.75501)
12 Redfish pass -82.200681 26.553486 12 -82.419777 26.430613 -81.981354 26.676300 POINT (-82.20068 26.55349)
13 Sao Fransisco -48.713194 -26.312026 11 -49.151459 -26.557822 -48.274612 -26.065419 POINT (-48.71319 -26.31203)
14 Senegal delta -16.667135 13.766087 11 -17.105713 13.499143 -16.228867 14.032679 POINT (-16.66713 13.76609)
15 Severn estuary -2.680262 51.595222 11 -3.118744 51.424474 -2.241898 51.765715 POINT (-2.68026 51.59522)
16 St Michel -1.510683 48.633302 13 -1.620312 48.587850 -1.401100 48.678607 POINT (-1.51068 48.63330)
17 Wadden sea 5.195870 53.300644 10 4.319000 52.970973 6.072693 53.627539 POINT (5.19587 53.30064)
18 Wax Lake -91.448385 29.498476 13 -91.557999 29.438701 -91.338787 29.558228 POINT (-91.44839 29.49848)
19 Willapa bay -124.071405 46.682020 12 -124.290733 46.587653 -123.852310 46.776082 POINT (-124.07140 46.68202)
20 Yangtze river 121.878402 31.289810 9 120.124512 30.346806 123.631897 32.224419 POINT (121.87840 31.28981)
21 Yellow river 119.173973 37.754205 11 118.735428 37.536955 119.612274 37.971267 POINT (119.17397 37.75420)
coastal_systems2.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook